%Matlab code for computing all ARMA(p,q) models ranging from p
%equal to 1 to 5 and q equal 0 to 5.  This code supports the manuscript "Putting the
%significance of spectral peaks on the level: implications for the
%1470-yr peak in Greenland δ18O" by P. Huybers, Feb 2022.  File were
%run on Matlab release R2021b and require the signal
%processing toolbox. 

clear;
addpath Toolbox;

%%load data
GISP2;
timeo=G(1:end,3)/1000;
Oo=G(:,2);
pl=find(timeo>=20 & timeo<=50);
O=Oo(pl);
time=timeo(pl);
mdt=mean(diff(time));
pl=find(timeo>=0 & timeo<=70);
Oo=Oo(pl);
timeo=timeo(pl);

%interpolate to uniform intervals
ti=time(1):mdt/10:time(end);
Oi=interp1(time,O,ti); 
Oi=smoothPH(Oi,11); 
time=ti(1:10:end); 
O=Oi(1:10:end);
dt=diff(time(1:2));

%%spectral analysis 
p_all=1:5;
q_all=0:5;
opts.nw=2;
dof=2*(2*opts.nw-1);
k=dof/2;
Vexp=psi(1,k);  %expected variance 

for ctp=1:length(p_all);
    for ctq=1:length(q_all); 
        opts.p=p_all(ctp);
        opts.q=q_all(ctq); 
        [P s coeffs]=specw(detrend(O),dt,opts.nw,opts);
        vP(ctp,ctq)=var(P.prewhiten); 
        F(ctp,ctq)=vP(ctp,ctq)./Vexp;    
        [ctp ctq F(ctp,ctq)],
    end;
end;
